Population extinction in a fluctuating environment 
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Environmental noise can cause an exponential reduction in the mean time to extinction (MTE) of 
an isolated population. We study this effect on an example of a stochastic birth-death process with 
rates modulated by a colored Gaussian noise. A path integral formulation yields a transparent way of 
evaluating the MTE and finding the optimal realization of the environmental noise that determines 
the most probable path to extinction. The population-size dependence of the MTE changes from 
exponential in the absence of the environmental noise to a power law for a short-correlated noise 
and to no dependence for long-correlated noise. We also establish the validity domains of the limits 
of white noise and adiabatic noise. 
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Extinction of species after maintaining a long-lived 
self-regulating population is a dramatic instance of a 
large fluctuation. Its origin is in the intrinsic discrete- 
ness ( "quantization" ) of individuals, and in the ran- 
dom nature of birth-death processes [H, Q ■ Understand- 
ably, the extinction phenomenon has attracted much at- 
tention from population biologists and epidemiologists 
3J. As birth-death processes are intrinsically far-from- 
equilibrium, they are also of much interest to physics 
[J [2} • Birth-death processes often occur in time- varying 
environments. Understanding the impact of environmen- 
tal noise on the mean time to extinction (MTE) of a 
species is both important [i| and interesting. Early mod- 
els assumed that the environmental noise, which modu- 
lates the birth and /or death rates of the species, is delta- 
correlated in time p, Q . More recently numerical simula- 
tions of the effects of a finite correlation time of the noise 
were performed by many population ccologists, see e.g. 
0- Not surprisingly, the simulation results provide only 
a partial understanding of the complex and rich inter- 
play between the nonlinear kinetics and intrinsic (demo- 
graphic) stochasticity of the population on the one side 
and the magnitude and spectral/correlation properties of 
the environmental noise on the other. 

In this Letter we formulate a theoretical framework for 
this problem by considering a prototypical example of 
a single-species stochastic birth-death process with rates 
modulated by a "red" noise: a positively correlated Gaus- 
sian noise with given magnitude and correlation time. We 
evaluate the MTE analytically and find that the qualita- 
tive and quantitative details of the exponential reduction 
of the MTE by the environmental noise are very sensi- 
tive to the noise color. It was discovered by Leigh |1, @] 
that white environmental noise changes the dependence 
of the MTE on the metastable population size from an 
exponential to a power-law with a large exponent. Here 
we show that a colored noise changes this exponent, re- 
ducing it at a fixed noise magnitude. For a very long 
correlation time of the environmental noise, where we 



develop an adiabatic theory, the MTE becomes indepen- 
dent of the population size for a strong enough noise. We 
also establish the validity domains of the limits of white 
noise and adiabatic noise. 

The distinct effect of the environmental noise on the 
MTE comes from special realizations of the noise which 
affect the birth and/or death rate in an optimal way. The 
optimization involves a statistical "cost" of a given varia- 
tion of the reaction rates along with a "gain" due to a fa- 
cilitated extinction. We find that the optimal realization 
of noise (ORN), which determines the least improbable 
path to extinction, changes considerably as the noise cor- 
relation time is varied. For a short-correlated noise the 
ORN has a form of a sudden "catastrophe" , bringing the 
reaction rates, for a certain period of time, to such val- 
ues that cannot sustain a steady-state population. For a 
long-correlated noise the ORN merely gradually reduces 
the population size. While not directly causing extinc- 
tion, it makes a fatal demographic fluctuation much less 
improbable. The ORNs in different intermediate regimes 
(depending on the rescaled noise magnitude and correla- 
tion time) can be found numerically. 

To be specific we consider a continuous-time birth- 
death process in a population of n species with the birth 
rate A„ and death rate [i n given by 

A„ = - (fj. + r - an) , //„ = - (// - r + an) . (1) 

For time-independent rate constants /i, r and a (we as- 
sume r < /j,), this is a symmetrized version of the lo- 
gistic Verhulst model: a well-studied model of popula- 
tion dynamics, see e.g. Q. The terms nonlinear in n 
describe, at a > 0, competition for resources which lim- 
its the exponential population growth. As a result, the 
rate equation n — rn — an 2 predicts, at r > 0, a sta- 
ble population of the average size n = K = r/a^>l 
which sets in after the relaxation time t r = 1/r. Demo- 
graphic noise, however, makes the "stable" population 
metastable. The population actually goes extinct, as a 
large fluctuation ultimately brings it to the absorbing 



state n = 0. Large fluctuations are rare and therefore 
statistically independent. As a result, the long-time sur- 
vival probability obeys Poisson's law 



J2P n (t) = l-P (t) = e- 



t/T 



(2) 



where P n (t) is the probability to find n individuals at 
time t, and r is the MTE. It is a well-known result (that 
we will reproduce shortly) that r scales exponentially 
with the average population size K, see e.g. Ref. [8|. 
In the limit of r <C fi, that we will be interested in, one 
obtains with exponential accuracy: 



To oc exp(r Kj /i) 



(3) 



where we have assumed rK / 'fj, ^> 1. 

Environmental noise manifests itself as a time- 
modulation of the birth and/or death rates. We will 
assume a modulation of the parameter r: 



r -> r(t) = r - £(t) 



(4) 



where £ (t) is a "red" (positively correlated) Gaussian ran- 
dom process with zero mean, variance v <C /i 2 and corre- 
lation time t c . For convenience, we choose the Ornstein- 
Uhlcnbcck noise defined by the correlator (£(i)£(i')) = 
■yg-lt-* l/ f c The statistical weight of a given realization 
of this noise is P[£{t)} oc exp{— 5[£(i)]}, where 



(5) 



The environmental noise does not change the Poisson 
character of the survival probability, Eq. @. Unless the 
noise is too weak, however, it exponentially reduces the 
MTE. We found that the MTE T£, reduced by the noise, 
can be expressed in terms of the unperturbed MTE r 
and two dimensionless parameters: the rescaled noise 
variance V — vK/(rfi) and the rescaled noise correlation 
time T — t c /t r = rt c : 



lnr £ = F(V, T) lnr 



(6) 



where the function F(V,T) describes various parameter 
regimes summarized in Fig.[TJ Importantly, each of these 
regimes is also characterized by an optimal realization of 
the noise (ORN) which causes population extinction with 
the highest probability. Not only the different regimes 
have exponentially different MTEs: they also feature 
qualitatively different ORNs. 

Our theory, which leads to Eq. Fig. [T] and other 
results, starts from the master equation for the time- 
dependent probability distribution function P n (t): 



Pn — Ki-lPn-l — (A„ + f^n)Pn + Mn+l-Pn + l 



(7) 



with the birth and death rates given by Eq. (TJJ . One can 
show that, for K 3> 1 and r <C n, this master equation 
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FIG. 1: Various regimes of extinction on the plane of rescaled 
parameters V and T. The dashed lines are schematic borders 
of the adiabatic, Eq. QB) , and white-noise, Eq. (|14|l limits. 
The dotted line is the border of the weak-noise, Eq. (|17[1 . 
regime. The shaded area is a crossover region. 



can be accurately approximated by the Fokkcr-Planck 
equation, derivable by a standard procedure of van Kam- 
pen system size expansion 0, 0] ■ Switching to the contin- 
uous notations n — * q, one can write the Fokker-Planck 
equation as P = HP, with the linear differential operator 



H{q,p) = -fq + p{rq - aq 2 ) 



(8) 



Here p = —d q so that [q,p] = 1. In the presence of envi- 
ronmental noise, see Eq. (j4|), one obtains the Hamiltonian 
H ( {q, P , t) = H(q,p) - £(t)pq. 

The evolution operator U(qf,tf,qi,ti) of the Fokker- 
Planck equation can be represented as a path integral 
over time-dependent trajectories q(t) and p(t). Below we 
discuss the boundary conditions for such trajectories for 
the case of population extinction. Eventually the evolu- 
tion operator must be averaged over realizations of the 
environmental noise, resulting in 



(U) = / ViVqVp e 



-Sl£]-fdt\pq-H(q,p)+Spq] 



(9) 



where p(t) and H(q,p) are understood as "classical" vari- 
ables, rather than the operators. 

Rare events in general and population extinction in 
particular are described by "classical" trajectories ac- 
cumulating a large action (and therefore having expo- 
nentially small probabilities). For this reason the cor- 
responding path integral can be evaluated using the 
saddle point approximation near the most probable (or 
rather least improbable) trajectory, describing a given 
rare event. Such an optimal trajectory is determined by 
the variation of the exponent in Eq. ([9]) over q(t), p(t) 
and £ (t) . The variation over £ yields the ORN which de- 
termines a given rare event with the highest probability. 
Executing this program, one arrives at the following set 
of classical equations of motion for q(t), p(t) and £(£): 
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dH 



dp 

t 2 J - £ = 2vt c pq . 
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The boundary conditions, corresponding to extinction 
of the metastable population of average size K, are 
q(t = -oo ) = K , q(t = +00) = 0, and f (t = ±00) = 0. 
The conditions for £ follow from the fact that the ORN 
must have a finite duration. Indeed, there is no need in 
environmental variations well before a large fluctuation 
starts and well after the population goes extinct. With 
exponential accuracy, the extinction probability of the 
large fluctuation is given by the full action, see Eq. J9]), 
calculated on the solution of Eqs. (fT0|) and 

In the absence of environmental noise, £ = 0, Eqs. (JTUJ) 
admit an integral of motion: H = const. Then it is easy 
to see that the trajectory obeying the proper boundary 
conditions has H = and is therefore implicitly given 
by the relation (fJ,/2)p + r — aq = 0, see Eq. (J5J and 
Fig. [2] Calculating the action along this trajectory one 
finds S = f^pdq = r 2 / (fia) — rK/fi which yields the 
extinction time ([3]). Solving for q(t) and p(i) for this 
trajectory, one finds the optimal path to extinction: 
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(12) 



In what follows we analyze, in various limits, the solu- 
tions of Eqs. (dU]) and (jlip in the presence of environ- 
mental noise. 

Short- correlated noise. Here the term t 2 £(t) in Eq. ([TTJ) 
can be neglected, and the ORN becomes enslaved to the 
dynamics of q and p: £(t) ~ —2vt c pq. As a result, 
Eqs. pop become Hamiltonian equations of motion with 
the effective Hamiltonian 



H v (q,p) = H{q,p) 



,22 
vt c p q . 



(13) 



[The same conclusion follows from the gaussian inte- 
gration over in Eq. ^ with the white-noise action 
S[£\ = J dt £ 2 / (Avtc) .] Now, H v is an integral of motion 
of Eqs. (fTU)) and, by virtue of the boundary conditions, it 
takes the value H v =0. As a result, the extinction pro- 
ceeds along the line (p/2)p + r — aq + vt c pq = 0, depicted 
in Fig. [5J Evaluating the action S = f^pdq along this 
line, one finally arrives at Eq. ([6]) with 



F{V,T) 



1 

VT 



1 + 2VT 
2VT 



ln(l + 2VT) - 1 



(14) 



As the (effectively) white noise is fully characterized by 
the product vt c , F only depends on the product VT. 

For a weak noise, VT < 1, Eq. (TJJ} yields F ~ 
1 - 2VT/3. The corresponding reduction of the MTE 
is still exponentially large as, according to Eqs. ([3]) and 



T e 



-2vt c rK 2 /3/T 



<C tq. However, the most dra- 



matic reduction of the MTE is predicted, in the spirit 
of the pure white-noise result 0, [B| , in the strong-noise 
limit, VT > 1. Here F ~ ln(VT)/(VT), and one obtains 



T 6 oc {vtcK/nY 1 ^ 



One can see that the exponential scaling of the MTE 
with the population size K, cf. Eq. ([3]), gives way here 
to a power law of K with a large exponent. To clearly 
see the origin of this qualitative change in the MTE, let 
us find the ORN leading to Eq. (Ti"5|) . The logarithmic 
term in Eq. (TTJJ comes from the hyperbolic part of the 
extinction trajectory, see Fig. [2l where vt c pq ~ — r. Here 
£(i) ~ —2vt c pq ~ 2r ~ const, and therefore q ~ — rq. 
The latter equation describes the population size decay 
from the initial value K down to u/(vt c ). At this scale 
[and at time t = t r ln(Kvt c / fx)], the demographic noise 
takes over the environmental one, see Eq. (fT3]) . As a 
result, the ORN of the short-correlated environmental 
noise is a catastrophic event [9|], where the parameter 
r > suddenly drops to — r and keeps this value for a 
logarithmically long time i » t r , see Fig. O The MTE 
([i"5j) merely reflects the statistical weight of this ORN. 
This argument also shows that the validity of Eq. (fT4|) 
requires a less restrictive condition then t c <C t r . Indeed, 
it suffices to demand that t c <C t = U ln(VT), see Fig. [T] 
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FIG. 2: Zero-energy trajectories of the Hamiltonian H (the 
dashed line), H v (the dotted line), and both H and H v (the 
solid lines). The shadowed area is the extinction action for 
the short-correlated environmental noise, leading to Eq. (|14|) . 



Long- correlated noise. Here an adiabatic theory can 
be developed. The rare fluctuation, causing extinction, 
takes time about t r , see Fig. [3] As the environmental 
noise changes on a much longer time scale t c , the ex- 
tinction fluctuation samples an almost constant value of 
the noise £(0) = £0, to be determined below. The ef- 
fective parameter r is therefore equal to r — £0 and is 
constant. Therefore, the corresponding extinction rate 
is ~ exp[— (r — £ ) 2 /(na)), cf. Eq. ([3]). Now we notice 
that the right hand side of Eq. (JTTJ) vanishes everywhere 
except in a small time window \t\ < t r <C t c . As a result, 
the solution of Eq. dHJ) for the ORN is £(i) ~ £oe" |t|/tc . 
Using it in Eq. (O, we find the statistical weight of the 
ORN to be ~ exp[— $/(2v)]. Finally we need to find 
the optimal value of £0 by optimizing the extinction rate 
against the statistical weight of the ORN. This is done by 
finding the minimum of £q/(2d) + (r — £o) 2 /(M a ) which is 
achieved at £0 = r[l + fia/ (2v)]^ 1 . The minimum action, 
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2 /([j,a + 2v), yields the logarithm of the extinction time, 

with 



which is therefore given by Eq. 

F(V,T) = (1 + 210 -1 



(16) 



Notice that, for a strong long-correlated noise, V 3> 1 and 
T» 1. one obtains lnr^ = r 2 /2v which is independent 
of the population size K. 

When does Eq. (fTB"]) apply? It turns out that, for a 
strong long-correlated noise, the condition T » 1 gives 
way to a more restrictive one. Indeed, when deriving 
Eq. (fl~6|) we assumed that r(t) ~ r — £oe~'*'/ tc does 
not change during the relaxation time t r . This requires 
r'(0)t r <C r(0) and leads to the condition T 3> max(l, V), 
shown in Fig. Q] as the border of the adiabatic regime. 

Weak noise. Here one can solve Eqs. ([TOjl and (jTTJ) 
perturbatively. This is equivalent to performing the in- 
tegration in Eq. over the unperturbed extinction tra- 
jectory, Eq. (|12p [10|| . The gaussian integration over the 
noise is done by going to the frequency space, and we 
obtain 



F(V,T) = 1-4V / — 



T 



2tt sinh 2 ttlo 1 + K) ; 



(17) 



For a short-correlated noise, T -C 1, this expression yields 
F = 1 - 2VT/3 in agreement with the limit of VT <C 1 
of Eq. (fT4"]) . On the other hand, in the adiabatic limit, 
T > 1, Eq. (JT7J) yields F = 1 - 2V in agreement with 
Eq. (flU)) at V -C 1. These arguments provide the border 
of the weak-noise result (TTT)) . depicted in Fig.[TJ We stress 
that even a relatively weak noise causes an exponentially 
large reduction of the MTE. 

Equation (| 1 7[) shows that, for a weak environmental 
noise, there is only one relevant scale for the noise cor- 
relation time: T ~ 1. The situation is more complicated 
for a strong noise, V 3> 1. As was discussed above, the 
adiabatic regime holds when T 3> V, whereas the ef- 
fectively white-noise regime holds for T <C lnT^. In the 
crossover regime, \nV < T < V (see Fig.Q}, the function 
F changes by a numerical factor of order unity. The MTE 
in the crossover regime can be inferred from a numerical 
solution of Eqs. (fT0|) and (| 1 1 1) which is quite straightfor- 
ward. 

To conclude, we have evaluated the reduction of the 
mean time to extinction (MTE) of an isolated popula- 
tion caused by environmental noise. We have also es- 
tablished the validity domains of the limiting cases of 
white and adiabatic noises. Even a relatively weak en- 
vironmental noise causes an exponential reduction of 
MTE. A strong noise brings about qualitative changes 
in the scaling of MTE with the metastable population 
size K . While MTE scales exponentially with K in the 
absence of environmental noise (or if the environmental 
noise is weak) , the scaling changes to a power law in the 
limit of a strong short-correlated noise, and becomes K— 
independent in the limit of a strong long-correlated noise. 




FIG. 3: Optimal realizations of the environmental noise in 
the limit of short (the dashed line) and long (the solid line) 
correlations of the noise. The duration of the "catastrophe" 
for the short-correlated noise is t ~ t r ln(Kvt c /fi). 



The optimal realization of the environmental noise, which 
results in the population extinction with the highest like- 
lihood, also differs qualitatively in these two limits. For 
a short-correlated noise the ONR has the form of a sharp 
"catastrophe" which, for a logarithmically long time, in- 
terchanges the birth and death rates of the system. For a 
long-correlated noise the ONR is a slow suppression of the 
birth rate down to a positive value. It is still debated in 
population biology "whether and under which conditions 
red noise increases or decreases extinction risk compared 
with uncorrelated (white) noise" , see Ref. 0(b). We 
hope that the analysis presented here will help resolve 
this and related issues. 

We thank M. Dykman for enlightening discussions. 
A. K. was supported by NSF grant DMR-0405212. B. M. 
was supported by the Israel Science Foundation (grant 
No. 408/08). 



[1] 
[2] 
[3] 



[5] 
[6] 
[7] 



[8] 
[9] 
[10] 



[11] 



C.W. Gardiner, Handbook of Stochastic Methods 
(Springer Verlag, Berlin, 2004). 

N.G. van Kampen, Stochastic Processes in Physics and 
Chemistry (North-Holland, Amsterdam, 2001). 
M.S. Bartlett, Stochastic Population Models in Ecoloyy 
and Epidemiology (Wiley, New York, 1961). 
S. R. Beissinger and D. R. McCullough (Editors), Pop- 
ulation Viability Analysis (University of Chicago Press, 
Chicago, 2002). 

E. G. Leigh, Jr. J. Theor. Biol. 90, 213 (1981). 

R. Lande, Amer. Naturalist 142, 911 (1993). 

(a) K. Johst and C. Wissel, Theor. Popul. Biology 52, 91 

(1997); (b) M. Schwager, K. Johst, and F. Jeltsch, Amer. 

Naturalist 167 (2006). 

I. Nasell, J. Theor. Biol. 211, 11 (2001). 

M. Assaf, A. Kamenev and B. Meerson. larXiv:0803. 04381 
A conceptually similar weak-noise theory has been re- 
cently developed by Dykman et al. in the context of the 
impact of random vaccination on decease extinction [ll| . 
M. Dykman, I. B. Schwartz, and A. S. Landsman, Phys. 
Rev. Lett. 101, 078101 (2008). 



